Adiabaticity in semiclassical nanoelectromechanical systems 
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We compare the semiclassical description of NEMS within and beyond the adiabatic approxi- 
mation. We consider a NEMS model which contains a single phonon (oscillator) mode linearly 
coupled to an electronic few-level system in contact with external particle reservoirs (leads). Us- 
ing Feynman- Vernon influence functional theory, we derive a Langevin equation for the oscillator 
trajectory that is non-perturbative in the system-leads coupling. A stationary electronic current 
through the system generates nontrivial dynamical behavior of the oscillator, even in the adiabatic 
regime. The 'backaction' of the oscillator onto the current is studied as well. For the two simplest 
cases of one and two coupled electronic levels, we discuss the differences between the adiabatic and 
the non-adiabatic regime of the oscillator dynamics. 

PACS numbers: 71.38.-k, 73. 21. La, 85.85.-l-j 



I. INTRODUCTION 

Nanoelectromechanical systems (NEMS) enable the 
detailed study of the interaction between electrons, tun- 
neling through a nano-scale device, and the degrees of 
freedom of a mechanical system. The electronic current 
affects the mechanical system and vice versa. The dimen- 
sions of systems used in recent experiments range down 
to scales, where the observation of fundamental quan- 
tum behavior for a comparatively macroscopic object is 
possibleii^. 

The influence of strong clcctron-phonon coupling in 
molecules or suspended quantum dots yields highly in- 
teresting effects^"—, like the Franck-Condon blockade, 
where the influence of the mechanical system sup- 
presses the electronic current^i^, or switching in molecu- 
lar junctions^. 

In non-equilibrium, a standard approach to solve the 
dynamics of NEMS is to do perturbation theory in a tun- 
nel Hamiltonian, which has been successfully done up to 
the co-tunneling regime^. This is a well-explored path, 
where the dynamics is described by master equations or 
generalizations thereof in Liouville space. Many interest- 
ing physical results can be obtained via this approach, 
e.g. avalanche- type molecular transport^ or laser-like 
instabilities^. These methods produce good results in the 
range of high bias, where non-Markovian effects can be 
neglected. To gain access to the small bias regime, one 
has to work perturbatively in the system-oscillator in- 
stead of the system-leads couplingiS. Alternatively, if the 
oscillator is treated in a semiclassical regime, Feynman- 
Vernon influence functional techniques are suitableii"— . 

In a previous work, we combined a semiclassical analy- 
sis with an adiabatic approacbi^, where we assumed the 
oscillator movement to be slow compared to the electrons 
which are jumping through the system. The condition 
r ^ wq then followed from this approximation, where T 
denotes the tunneling rate of the electrons and wo is the 
oscillator frequency. The most interesting results, such as 
negative damping and limit cycles, were achieved in the 
regime where the oscillator and the electrons act on the 



same time scale. But in the latter regime, the adiabatic 
approach was at the limits of its validity. 

In recent publications different approaches, e.g. based 
on scattering theorj*iiii^, were used to go beyond the adi- 
abatic approximation or to verify the range of validity 
for this approachi^^. In this paper, we go one step fur- 
ther and present numerical results for a completely non- 
adiabatic approach. We apply our method to two sim- 
ple NEMS models, the single and the two-level system. 
We work in the semiclassical regime where an expansion 
around the classical path is performed. The advantage 
of this method is that we are non-perturbative in the 
system-leads coupling, because the exact electronic solu- 
tions are included. In the non-adiabatic approach, we can 
treat the oscillator and the electrons on the same time 
scale without further constrains. Therefore, we modify 
our adiabatic path-integral approach, obtaining an ex- 
plicit time-dependent perturbation. As a consequence, 
wc have to calculate system quantities numerically in a 
full time dependent manner. This allows us to critically 
assess the validity of the adiabatic results. By compar- 
ing the outcomes of the adiabatic and the non-adiabatic 
approach we find a qualitative accordance. The same 
features arise in both approaches and the stationary re- 
sults predominantly coincide, but there is no quantitative 
accordance in the time-dependent regimes. The differ- 
ences increase together with the complexity of the fo- 
cused NEMS model. 

This paper is organized as follows. In Sec. |II] we in- 
troduce the general model, followed by the derivation of 
the stochastic equation of motion for the oscillator dy- 
namics. In Sec. mil we present the results for a single 
resonant level system, comparing the adiabatic and the 
non-adiabatic regime. In addition to the phase space 
trajectories, we calculate the resulting electronic current 
through the system. In Sec. IIVI we consider a two-level 
system, where the dynamical behavior of the oscillator 
exhibits non-trivial effects. Thereby we show that the 
oscillations of the mechanical subsystem lead to an oscil- 
lating current with fixed frequency. 



II. MODEL 

Our total Haniiltonian is a sum of an electronic system 
Tic, a single oscillator with a spatial degree of freedom 
g in a harmonic potential Hose and a linear coupling be- 
tween the oscillator and the electronic system 



n — ric + rio 



Fq, 



(1) 



in which F denotes an electronic force operator. The 
electronic system itself consists of a few electronic levels 
which are connected to two macroscopic leads. The latter 
are considered as two Fermi seas with chemical potential 
A*aGL,R and temperature T. Furthermore, the electronic 
part provides a non-equilibrium environment for the sys- 
tem. In our formalism we work non-perturbatively in the 
system-leads coupling, assuming arbitrary coupling and 
a finite bias regime without constrains. 

The single oscillator with momentum p, position q and 
mass m is described by a parabolic potential 



rto 



1 

2m 



P 



-muj^f, 



(2) 



whereby loq equals the oscillator frequency. The oscillator 
potential will be modified by the electronic environment 
and will exhibit multi-stabilities due to the electronic 
forcesi^. In this paper, the reduced Planck constant is 
set to unity {h ~ 1). 



A. Influence functional 

We want to focus on the oscillator's dynamics, which is 
described by the reduced density matrix posc (t) , obtained 
from the total density matrix x(t) by tracing out the bath 
degrees of freedom, posc(t) ~ trBx(0- The propagation 
of the reduced density matrix in time can be written as a 
double path integral over qt (forward) and q'^ (backward) 
weighted by the Feyman- Vernon influence functional^ii^ 



T[qt;q't]=tTB{U^q't]U[qt]), 
containing the time evolution operator 



U[qt 



-i fdt' \H^-Fq,, 
Te » 



(3) 



(4) 



The time-ordering operator T arranges operators with 
later times to the left. 

We want to perform an expansion around the classi- 
cal path. Therefore we transform to center-of-mass and 
relative path variables 



qt 



1 
xt + -jVt, 



1 



% 



Xt 



:yt, 



(5) 



where the variable yt can be interpreted as the quantum 
fluctuations around the classical path. We introduce an 



interaction picture 

U[qt] 



Te " 



^ _ i f dt' iF(t')y,/ 

= U[xt]U[yt]; U[yt]^Tei ,(6) 

where the term with the off-diagonal path yt is regarded 
as a perturbation and F(t) = W [xt] F U \xt\- Inserting 
Eq. ^ into the influence functional, Eq. ([3|), leads to 

F[quq[] = trB{U^[~yt]U^[xt]U[xt]U[yt]) 

^ tMUH~Vt]U[yt]). (7) 

Expanding this term to second order and performing a 
cluster expansion^, we finally obtain 



with the influence phase 



-$[a;t;i/t] 



(8) 



* [xuVt] = - i / dt'f{t')yt, + dt' dsC{t',s)yt'ys, 



(9) 
and the force correlation function 

C(t', .) = trB { [m fit')) [f{s) /(.)) } 

^{5F{t')6F{s)). (10) 

The force term f{t) = {F{t)) depends on the center of 
mass path xt. The influence phase, Eq. ©, can be re- 
garded as a cumulant generating functional for the force 
operator correlation functions. The term of quadratic or- 
der in yt describes the Gaussian fluctuations around the 
classical trajectory which is determined self-consistently 
in our approach. Higher order terms in yt (corresponding 
to non-Gaussian noise) describe higher quantum fluctu- 
ations which are neglected here. 

B. Langevin equation 

To second order in yt the double path integral for the 
reduced density matrix describes a classical stochastic 
process for the diagonal path xt that is defined by a (non- 
adiabatic) Langevin equation 



mxt + V^^^ixt) - f [xt] = 6 



(11) 



with V^g^{xt) ~ mujQXt and a Gaussian stochastic force 
^t that has a correlation function {^t'£,s) = C{t',s). 
Eq. pT|) is the starting point for our non-adiabatic cal- 
culations. Note that the force / [xt] is a complicated 
functional that contains the full time-dependence of the 
position operator xt- 

In the adiabatic approximation, a Taylor expansion for 
the center of mass variable is performed (xj w xo + t*xo), 
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FIG. 1: Phase space portraits resulting from the (non-) adi- 
abatic approach in units of luq and lo with the parameters 
r — 1.4^0, £d ~ S.Oluq and g = 2.45 at zero tempera- 
ture. Row A depicts the adiabatic resuhs, additionally the 
case without friction is plotted (dotted lines). Row (B) de- 
picts the non-adiabatic results. The bias voltage (symmet- 
ric choice) is increased from left to right, explicit values are 
Vbias/wo = 0.5/2.5/5.0. 



leading to an interaction picture with respect to Hq = 
He — Fxq and a perturbation y[g](t) = —F{txQ + ^yt)- 
Consequently, the expectation value of the force operator 
in the adiabatic interaction picture can be calculated for 
fixed xq . Additionally, due to the second term iig of the 
Taylor expansion, an explicit friction term arises in the 
(adiabatic) Langevin equation, 

mxt + V^,,{xt) + XtD [xt] - J\xt\ = 6- (12) 

The friction term D [xt\ can be interpreted as the first 
adiabatic correction term. 

In contrast, in the non-adiabatic Eq. (jlip . all higher 
orders of the Taylor expansion are included and the first 
challenge is to calculate the force term / [xt] considering 
the full time-dependence of xt- In all our calculations 
presented in this paper we neglect the stochastic fluctu- 
ations (^t =0). 



III. SINGLE RESONANT LEVEL 

The Hamiltonian for the single resonant level 
(Andcrson-Holstcin model - AHM) is 



^-E 



^ka^l^Q^^k 



E (^^"S: 



d + h.c. 



■ d^d+Uo 



ka 



ka 



(13) 

with the abbreviation e ~ Sd — Xq, containing the energy 
of the local level Sd and the linear coupling to the os- 
cillator. Thereby, A equates the coupling strength. The 
operators d/d^ correspond to the dot, and the lead oper- 
ators Cka/Cka annihilate/create an electron in the a-lead 
with energy Ska and momentum k, whereby a = L, R de- 
note the left/right lead. Transitions of electrons between 



the dot and the a-lead are possible with the amplitude 
Vka- Here, the force operator becomes F = Xd)d = Xhd 
and the self-consistent equation of motion for the expec- 
tation value of the oscillator coordinate Xt (neglecting 
stochastic fluctuations) reads 



TOXt + V^^^ixt) - \N [xt] =0 
N[xt] = tvB{U^[xt]hdU[xt]). 



(14) 



The occupation number of the local level N [xt] is calcu- 
lated with the lesser Green function 

N[xt] ^{d\t)d{t)) = -iG<{t,t), (15) 

which is obtained fronr the Keldysh equation 

G<{t,t) = fdh f dt2 G'^itM) ^^{hM) G"(t2,i), 



containing the advanced/retarded Green function 

-^ijdt"[e{t")^t^] 

G"-'"(i,t') = Tie(±fTi) e *' 



(16) 



(17) 



with e{t) = Ed—Xxt- Thereby, 8 designates the Heaviside 
step function. In the time dependent case the lesser self 
energy^! reads 

S<(fi,i2) = z^ /^e— (*^-*=)/„(c.)r", (18) 



whereby /^(w) denotes the Fermi function. We assume 
constant tunneling rates r„ = 27r^^ |VfcQp(5(u; — eua) = 
r/2, i.e. the left and the right tunneling rate are equal 
with r = Fl + Fr. As a result, for the time dependent 
occupation we obtain 



ivH= Y. T^ I ^ Uu^) \A{uM' 



(19) 



aeL,R 

with the spectral function 



-ifdt" (e(t")-w-iE) 

A{uj,t) = -i / dt' e *' 



(20) 



to 



To solve the equation of motioir without airy further ap- 
proximation and expansion, we transform Eq. (|20p into 
the differential equation 



A{uj,t) ^-i-i [ed-Xxt-uj-i- ] A{uj,t). (21) 



For the numerical integration a trapezoidal rule for dis- 
crete functions is applied. Thus we solve Eq. ([?!]) to- 



gcthcr with the system 



Xt =—pt 
m 



Pt 



^osc+aj: i^A, 



QeL,R 



N-l 



Y^ [/cK+l) |^K+l,t)P+/„K) \A{uJn,t)\' 



n=0 



(22) 



with Aw = \ujn — wo|/iV, in which N equals the number 
of points of the discretization scheme. 

Figurc[T]depicts the results in the oscillator phase space 
for different bias values. All data are obtained for a small 
tunneling rate F = IAljq, close to the limit of validity of 
the adiabatic approach which requires T ^ loq. In or- 
der to obtain correct physical units we introduce the di- 
mcnsionless coupling parameter g = A/(7Tia;gZo), whereby 
^0 = l/y/muiQ equals the oscillator length. Row A shows 
the adiabatic results. The dotted lines correspond to 
the adiabatic case without the first adiabatic correction 
term. Here, the trajectories run about the fixed points of 
the system. By varying the applied bias, the number of 
fixed points changes and in the case of high bias only one 
fixed point survives. Turning on the friction (first adia- 
batic correction term) leads to the solid line results in the 
graphs of row A. Here, the centers turn into stable spi- 
rals and the trajectories end up in the fixed points. (For 
further explanation see^^). In row B the non-adiabatic 
results are depicted. After long times t we observe little 
variation to the adiabatic result, nevertheless all trajec- 
tories end up in the same fixed points. By comparing 
both approaches the largest differences emerge for small 
times. This qualitative good accordance can also be ob- 
served in the results for the electronic current. In the 
non-adiabatic approach, the current is obtained frontSi 



I„(t) 



-eF„ 



N[xt] 



did 



U{L0)llIi[A{LU,t)] 



(23) 



The imaginary part of the spectral function is negative 
and describes the current flowing from the left lead into 
the dot. While keeping the oscillator center of mass co- 
ordinate Xt = X fixed in Eq. (j23p when calculating the 
spectral function, the adiabatic current result is repro- 
duced. Starting from Eq. ((20)) we obtain 



A 



adiabatic 



(w,i) 



-i(ed->'2;-w-?:£)(t-to) 



1 



(Ed - \x 



^i) 



(24) 



For large times t — to, the first exponential can be ne- 
glected and the spectral function becomes stationary. 
Therefore the adiabatic current reads (zero temperature) 



It = e- 



47r 



arctan 



2(m: 



L — £d 



Xx 



— arctan 



r 

2(mr - £d + A.T 



r 



= -Ir. (25) 




FIG. 2: LEFT: Current for t — ^ oo as a function of Vbias 
for two different coupling parameters g and with F — lAcjo, 
£d = S.OtJo at zero temperature. The black symbols depict 
the non-adiabatic and the indigo solid line the adiabatic re- 
sults. For comparison the infinite bias result (dashed line) and 
the current without coupling (dashed-dotted line) are plot- 
ted. The dotted lines mark the hysteresis like regimes, where 
two clearly distinct current channels exists. RIGHT: The up- 
per graph depicts the position dependent current Il(2;) for 
Vbias = 2.0a;o, the three fixed points are marked as crosses. 
The lower graph shows the left current for small times and 
Vbias = 1.0a;o. Remaining parameters are equal to the left 
graph and g — 2.45. 



This is a well known result and for the infinite bias case 
we obtain X]^ji = ±eF/4 as expected. Note, that in the 
adiabatic case the values for left and right current only 
differ in their sign. 

In the left graph of Figure [2] the stationary left current 
Il(^ -^ oo) is depicted for increasing bias and for two dif- 
ferent coupling parameters g. As explained above, the os- 
cillator trajectories end up in fixed points for large times. 
Hence, the current Il(^ -^ oo) becomes stationary and 
its value corresponds to a single level which is shifted by 
~gx* /Iq, whereby x* corresponds to a fixed point. Be- 
cause the system owns multiple fixed points, we obtain 
several current channels depending on the initial condi- 
tion. The dashed-dotted line in Figure [5] corresponds 
to the case without coupling to the oscillator {g = 0). 
For snrall bias, the current is small compared to the in- 
finite bias case (dashed line). There, the effective level 
e = Ed — gx* /lo is situated outside the transport window. 
The latter is also valid for the case without coupling, due 
to e ^ Sd ~ S.OujQ. The bias range for the current sup- 
pression is larger in the case of stronger coupling to the 
oscillator {g = 3.5). We obtain a hysteresis like shape for 
the current evolution, which is due to the multi-stability 
of the system. The coupling between the electronic and 
the mechanical system leads to a modified oscillator po- 
tential with additional minima. Switching between these 
states is possible and was theoretical proposed and stud- 
ied by several author o^'^i^^'^^ . 

The beginning and the ending of the hysteresis regime, 
where two current channels exist, are denoted by vertical 



dotted lines in Figure [2l For the non-adiabatic case the 
latter regime, where two current channels exist, differs a 
bit from the adiabatic case. 

The upper right graph of Figure [2] shows the current 
Ih{x) for the bias value Vbias = 2.0wo- Here, three fixed 
points occur, denoted by a cross. For (xj) ~ g/2lo the 
effective level is situated in the middle of the transport 
window and following from that the current is maximal. 
For the two other fixed points the effective level is again 
situated outside the transport window and the current is 
small. 

By comparing the adiabatic and the non-adiabatic sta- 
tionary currents, we can conclude that in the long-time 
limit only small differences exist. The differences are at 
their maximum for small times, which is clearly visible 
in the lower right graph of Figure [21 which depicts Zl (t) . 
There, the oscillations in the non-adiabatic case are much 
larger. 

The left and right time dependent currents differ for 
small times t in the non-adiabatic case. When we inte- 
grate II and Ir over all times the results coincide, so 
that there is no violation of current conservation. This is 
comparable to a periodically driven system with time de- 
pendent tunneling rates^l. The spectral function defined 
in Eq. (PO)) . is sensitive to small time differences t — to- 
For larger times the oscillator settles into one of the fixed 
points, whereas the spectral function becomes stationary 
and hence also the current. 



IV. TWO-LEVEL SYSTEM 

The model we are treating in this section consists of 
two single dot levels which are coupled by a tunnel bar- 
rier. Again we assume a coupling to a single bosonic 
mode. The total Hamiltonian is composed of the oscilla- 
tor part "Hose, cf. Eq. ([2]), the electronic part He and 
an interaction part which describes the coupling between 
the oscillator and the two dots. In contrast to the AHM, 
here the oscillator couples to the difference of the occu- 
pation numbers with the coupling strength A. The total 
Hamiltonian therefore reads 

n = n, + Hose - Xqidldi^ - 44), (26) 

containing the electronic part 



ka 



ka 



+ J2 ^o.dld^ + T, 4 dR + T; 4 dL, (27) 

a 

where h'aeh,R denotes the left and right dot energy lev- 
els. Here, Tc denotes the tunnel coupling matrix element 
between the two dots. Again, we obtain a Langevin equa- 
tion, Eq. pT|) . with the force term 



/H = A(a,)(i)^A[(nL)(i)-(nR)(t)], 



(28) 




FIG. 3: Phase space portraits for various tunnel couplings 
Tc, increasing from left to right. Upper row A: adiabatic 
results. Lower row B: non-adiabatic results. In graphs A2 
and B3 limit cycles (LC) appear. Explicit parameters are 
|Tcp = 0.2; 0.4; 1.0; 4.0 ljq. With the rate T = 2.0luo and for 
the chemical potentials we assumed fj.L = ^0 and /ir — — 5cjo. 
The dimensionless coupling constant is chosen as g = 2.5, and 
the internal bias voltage as V^nt = 5a;o, whereas i/l = —vn = 
eT/int/2. 



A. Time dependent occupation 

Calculating the time dependent occupations for the 
two level system using Green's functions is a challenge 
due to the complex dependencies and couplings of the 
systems operators. We choose a more direct way by us- 
ing the equations of motion technique, leading to a large 
system of coupled differential equations which have to be 
solved numerically. 

The Heisenberg equations of motion for operators of 
the dots and the leads {i>a{t) — i^a T Axt) yield 



Lit) = - I ( i>L(t) - ij) kit) - iT^Mt) + J2CkLit), 



where (ri,L/R)(i) implicitly depend on the oscillator coor- 
dinate Xt, cf. below. 



Mt) =-i(Mt)~ij)dRit)~ iTAit) + J2 CkRit), 



(29) 



where Ckait) = — iV^*Q,e~*'^'"°*Cfca(0) and the tunneling 
rate equals Fq = 2'Kj2k\^ka\'^Siu! — Ska)- Again, we 
assume constant tunneling rates Fl — Fr = F/2. 

The equations ([^ already include the solution for the 
inhomogeneous differential equation for the lead operator 
Cka ■ The tilde denotes the interaction picture introduced 
above, cf. Eq. ©. Hence, the effective time dependent 
energy level Vait) contains only the classical variable Xt- 

The differential equations for the corresponding dot 
annihilation operators are derived in a similar manner. 
Finally, one obtains an inhomogeneous system of coupled 
differential equations with time dependent coefficients. 
Multiplication with (5(aj — Ek.a) and summing over all k 



states leads to 



(^^) W = - 2 (^.)(^) + 2 Re 2 DRLit) 



+ / da;BLL(w,0- / duj'BYiR{uj',t) 



DRL{t) = i (.Mt) - Mt) + «0 DRL{t) - \Tcf (CT.)(i) 
+ dLo'BnL{uj',t)- I dujBl^{uj,t) 



Ba,ai^,t) = - i { Va{t) - U - i— ] Ba.aiuJ,t) 

r 

Ba,l3{uj,t) + -—fa{uj) 

47r 






(30) 



with the definitions: 



BnL{uJ,t) 



-T, VkK S{lo - ekn)e'''"'\dLiO)Mt)) 



Slr(w,0 = -T: VkL S{uj - ekL)e''''^H£l^{0)Mt)) 






i^RL(w,i) = *T, (4(i)dLW>- 



(31) 



The system Eq. ([50)1 is solved numerically together with 
the equations of motion for the expectation values for 
position and momentum operator 



1 

xt =—pt 
m 



ft=-K.c + A(a,)(t). 



(32) 



Hence, the phase space trajectories are obtained. 

In Figure [3] results for F = 2^0 are plotted. The upper 
row depicts the result for the adiabatic case including the 
first correction term. This so-called intrinsic friction term 
D [xt], cf. Eq. (fT2|). results from the non-equilibrium elec- 
tronic environment. The tunnel coupling |Tcp increases 
from left to right. Three fixed points appear in the range 
of small tunnel coupling (Al). There, the trajectories 
form stable spirals and run into the fixed points. For 
|Tcp = OAloq a limit cycle appears in the middle (A2). 
By further increasing \Tc\'^ the limit cycle turns into an 
unstable spiral(A3) and in the end only the left fixed 
points survives (A4). 

The second row (B) shows the results for the non- 
adiabatic approach. Qualitatively the same features 
emerge, as the appearance of the limit cycle and the 
stable spirals. Comparing the adiabatic and the non- 
adiabatic approach, we obtain quantitative differences, 
like the change of the middle fixed point into a limit 
cycle which happens at higher values of |Tcp as in the 
adiabatic case. The results differ most for small times, 
similar to the single level case (Scc lIII[) . By further de- 
creasing the tunneling rate F the differences between the 
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FIG. 4: LEFT: Current for i ^- cxd as a function of tunnel 
coupling iTcp. Al displays the results for XT^^ < 0.5luo and 
A2 the results in the range of 0.5 < \Tc\'^ /loq < 2.0. The 
symbols denote the non-adiabatic results. The diamonds cor- 
respond to the stable spiral situations, there the oscillation 
of the dynamical system disappears in the long-time limit 
and the current becomes stationary. Circles denote averaged 
current values for the limit cycle case, when the system per- 
forms periodic oscillations. The indigo sohd (dotted) lines 
depict the adiabatic results for the stable spiral (limit cycle) 
case. The dashed-dotted line depicts the current without cou- 
pling. RIGHT: The graph Bl shows the phase space results 
for \Tcf — 0.45a;o, here the radius for the adiabatic limit 
cycle is much larger than in the non-adiabatic case. Below, 
graph B2, depicts the corresponding time dependent left cur- 
rent, which oscillates as well. Explicit parameters are F = a;o, 
/iL = tJo and ^R = —5ljo- The dimensionless coupling con- 
stant is chosen as g = 2.5, and the internal bias voltage as 
Vint = 5uJo, whereas i^l = — fR = eV^nt/2. 



approaches increase. Results for a smaller tunneling rate 
will be presented in the next section. Sec. IIVBI there we 
calculate the electronic current for T ^ ojq. 

We also mention that the appearance of the limit cycle 
in this system is possible due to energy transfer processes 
between the electrons and the oscillator. In the stable 
spiral case, the influence of the electrons leads to damping 
of the oscillator. For example for i/l < %,, the electrons 
need energy to pass through the system. 

The requirements for a damped dynamical system to 
exhibit limit cycles is the additional appearance of pos- 
itive friction. For our system this means, that energy 
transfer processes occur which lead to the acceleration of 
the oscillator. The occurrence of positive friction is pos- 
sible when i>L > t'R and the electrons can transfer energy 
to the oscillator, cf J^. 



B. Current 

The current through lead a is derived via the Heisen- 
berg equations of motion and yields 



I„(t) 



^aNa [Xt] 



dujRc [Baail^it)] 



(33) 



The left graphs of Figure 2] depict the stationary current 
Ii^{t — > oo) as a function of the tunnel coupling |Tcp. In 
the adiabatic case and for small values of \Tc\'^ (Al) we 
observe two fixed points and one limit cycle leading to 
a tri-stable current. In the limit cycle case the current 
oscillates in time. The corresponding averaged current 
(dotted line) is not completely shown in (Al), due to 
the large values, since the current increases further until 
jTcP = 0.48ci;g. There, the limit cycle disappears and 
two fixed points remain until |Tcp = I.Swq (A2). We 
also obtain two fixed points in the small range of jTcp < 
0.03ojq, which is not dissolved in graph Al of Figure |31 

The current corresponding to the fixed point x* ~ 2/^o 
(solid line below g = case) increases approximately in 
the same fashion as in the case without coupling. In 
this regime the left effective level lays inside the trans- 
port window (t'L^R ~ =f2.5wo)- For large tunnel coupling 
one fixed point persists, a;* ~ — 2.3/Zo, and the corre- 
sponding current (lowest solid line) is strongly suppressed 
compared to the case without coupling. There, both ef- 
fective levels t'L,R — ±7.5a;o are clearly situated outside 
the transport window. Therefore tunneling through the 
two level system is rarely possible. 

In the left graphs of Figure S] the symbols denote the 
non-adiabatic current in the long-time limit. Here, the 
system has also two fixed points, but the limit cycle range 
is much larger, 0.35 < iTclV^o < 6-5. For IT^P < 0.35a;g 
wc observe two stable fixed points and by increasing the 
tunnel coupling the middle stable spiral turns into a limit 
cycle and the mechanical system performs periodic oscil- 
lations. For the latter case, the circles in Figure H] denote 
the averaged current. 

The non-adiabatic current corresponding to the mid- 
dle fixed point/limit cycle follows the result without cou- 
pling. As long as the fixed point x* ~ 0.03/Zo is stable the 
resulting effective level is approximately ?l.r — J^l.r as 
in the case without couphng. In graph Bl of Figure|31 the 
phase space trajectories are plotted for the case when the 
system performs periodic oscillations. The limit cycle, 
corresponding to the non-adiabatic results, runs in small 
cycles about the origin and following from that, the av- 
eraged current is similar to the case without coupling. In 
the adiabatic case the radius is much larger and the shape 
of the limit cycle is not smoothly circular. Hence, the cur- 
rent is much larger then in the non-adiabatic case (Al). 
In graph B2 the related time dependent current is de- 
picted. Frequency and amplitude differ strongly in both 
cases. The frequency of the current oscillations is equal 
to the oscillator frequency (non-adiabatic: w ~ 0.86ajo)- 

The current reaches its maximum when the distance 
between the left (right) effective level and the left (right) 
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FIG. 5: Time evolution for the current and the correspond- 
ing correlation functions for \Tc\ = uio in the non-adiabatic 
limit cycle case. The vertical dotted lines correspond to the 
maxima/minima of the effective level Ua = i^a T gxt/lo- The 
first row depicts the position for the effective left (Al) and 
right (Bl) level and the dotted line corresponds to the chem- 
ical potentials /iL ~ tJo and fin = —5ujo- The results for the 
left/right current are plotted in graph A2/B2, together with 
the occupation for the left/right level and the real part of 
the correlation function -Bll/rr- The time derivative of the 
level occupation is depicted in row three. Row 4 shows the 
result for the real (A4) and the imaginary (B4) part of the 
dot-dot-correlation function Drl. The real parts of -Blr/rl 
are depicted in row 5. The dimensionless coupling constant is 
chosen as g — 2.5, and the internal bias voltage as 'Knt ~ 5ijJq, 
whereas i^l ~ —I'n = el^nt/2. 



chemical potential is minimal (maximal). This is clearly 
visible in Figure [SJ where the time evolution for current 
and the different correlation functions defined above, cf. 
Eq. pT|) . are depicted. The first row shows the behavior 
of the oscillating effective levels, the value of the chemical 
potential is also plotted in these graphs (dotted line). In 
the first graph of the second row the current for the left 
lead is maximal when the effective level is minimal as 
mentioned above. The current decreases when the left 
level increases its distance to the transport window. 

The additional current peak (arrows in A2 of Fig. 
[5]) near the maximum of the left effective level does 
not appear in an adiabatic approach, where the cur- 
rent follows the position of the levels. This peak is 
related to internal coherent electronic oscillations be- 
tween the two dots. These oscillations are visible in 
the real part of the D^l function depicted in graph 
A4 of Figure [5] (denoted by dashed boxes), with fre- 
quencies that match the time dependent Rabi frequency 
u^Rit) = Vi^L(t)-%(t)+4|T,|2. 

In the adiabatic case, the time-resolved current for the 
right lead is equal to the current through the left lead 



with opposite sign. For the non-adiabatic case, right and 
left time- resolved currents are different, but their time- 
averages coincide. If we are in the long-time limit, and 
the system performs no oscillations, left and right current 
are equal. In contrast, in the limit cycle case we obtain 
a driven system leading to currents Il,r whose time de- 
pendence differ, since charge temporarily accumulates in 
the dots. 



V. CONCLUSION 

By comparing the adiabatic and non-adiabatic results 
for the single-level system we obtain a good qualitative 
agreement. In principle, the same features arise, as bista- 
bility and a hysteresis-like T — V characteristic are ob- 



served in both cases. The largest deviations are observed 
for small times, but in the long time limit the results 
predominantly coincide. 

For the two-level case the differences are much larger. 
Qualitatively we observe similar properties, but the 
quantitative predictions of the adiabatic approach do not 
match the results for the non-adiabatic system where the 
oscillator and the electrons act on the same timescale. 
The electron-oscillator interaction leads to multiple cur- 
rent channels like in the single-level system. Additionally, 
we observe limit cycles of the dynamical system leading 
to periodic oscillations of the current. In this regime, the 
system acts as a DC-AC-transformer. 
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